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1 .0 Introduction 


Considerable care Is taken In the design and construction of wings to 
ensure that the shape provides the required combination of lift and drag 
over the flight cycle and that the surface Is aerodynamical ly smooth. The 
presence of rain. Insect deposits or Ice can change the shape of the wing 
and Its surface finish and this paper examines the magnitude of the effects 
on lift and drag and describes the status of calculation methods which can 
provide a basic tool for their prediction. 

The problems associated with flying airplanes through heavy rain Include 
those associated with aerodynamic performance. It Is difficult to quantify 
these effects from flight experience since they occur usually together with 
other effects such as wind shear and downdraft. It Is known, however, that 
heavy rain can Increase the effective thickness of a wing and cause rough- 
ness which stems from drop Impingement and from waviness of the liquid 
film. These effects can, In turn. Influence the transition from laminar 
to turbulent flow and Increase drag while decreasing lift. They are sig- 
nificant at all angles of attack and can be Important at the higher angles 
associated with landing configurations. More will be said of this topic 
later In the course by Dunham. 

Knowledge of Insect contamination Is less, mainly because the likely conse- 
quences are small. It Is assumed that the contamination acts as distrib- 
uted roughness with maximum height In the region of the leading edge. The 
contamination tends to be removed at high speeds or In the presence of 
rain. 

The formation of Ice Is usually confined to the leading-edge region and 
again has maximum Importance at the high angles of attack associated with 
landing and takeoff. Deicing Is carried out where possible and can avoid 
or reduce the problem but, as Is known from recent accidents. Ice can form 
rapidly on the leading edge of wings and Intakes with considerable conse- 
quences for aerodynamic properties. Ice formation can also occur at cruise 
conditions and to an extent that lift Is reduced by an Important amount. 
It can be considered In two ways, the first where the effective shape Is 
changed and the second where the Ice acts as an equivalent sand-grain 
roughness although both can be Important In many circumstances. 
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I he following section examines the experimental evidence lor the effects 
of rain. Insects and Ice on airfoil performance and considers the extent 
to which the available Information can be Incorporated In a calculation 
method In terms of change of shape and surface roughness. It Is easy to 

envisage that a major shape change will have effects which can be described 
by the same procedures which led to the arrangement of the original air 
foil. In a similar manner, roughness can be Incorporated In the solution 
of boundary-layer or the Navler-Stokes equations provided the character- 
istics of the roughness are known. Thus the experimental knowledge of 
rain. Insects and Ice must be presented In the form of equivalent rough- 
ness. It Is also known that the environmental effects can affect the onset 
of transition and, since this can be Important with low Reynolds-number 
airfoils and with attempts to ensure laminar flow airfoils, this evidence 
Is also examined. 

The experimental Information has been used traditionally In the form of 
correlation equations and these are reviewed In Section 3. The advantage 
of these correlations Is that they can provide accurate representation 
within the limited range of the data but they are restricted by their lack 
of a physical basis for the equation. 

The fourth section of the paper considers the components of a method, based 
on more fundamental equations, to calculate the performance of airfoils as 
a function of shape, angle of attack and Reynolds number. One procedure 
Is described In greater detail and the ways of accommodating changes to the 
airfoil shape and surface roughness are considered. It Involves the numer- 
ical solution of conservation equations In differential form and has been 
used to obtain results which are presented In Section 5 and allow appraisal 
of the numerical features of the calculation method and of the extent to 
which It can predict the known environmental effects. 

2.0 Experimental Evidence 

The effects of rain have been examined In the wind-tunnel tests of Refer- 
ence 1 and more recently In References 2 and 3. The magnitude of the rain 
falls considered stem from arguments similar to those of Haines and Luers 
[4] who examined the records of the U.S. weather stations and concluded 
that the yearly mean maximum rainfall rate over a 60 sec period In the 
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eastern United States ranged from 150 to 250 mm/h. It Is expected that 
shorter term averages will achieve larger values and It should be noted 
that the record rainfall rate Is 1830 mm/h. Although these torrential 
rainfalls are uncommon. It Is desirable to know their likely consequences. 

The average thickness of the film of water on a 10m chord airfoil and fus- 
elage at zero angle of attack was calculated In Reference 4 and Is shown 
on lable 1. It Is unlikely that the distribution of film thickness would 
be uniform and Increasing angle of attack Is likely to lead to Increased 
thickness In the tral 1 Ing-edge region since the drag force between the air 
and water will decrease from around midchord. Thus, the tral 1 Ing-edge 
region can be expected to support film thicknesses considerably greater 
than those of lable 1, so that the effective shape of the airfoil can be 
altered by the rain to Imply an adverse pressure gradient In the aft part 
of the upper side of the airfoil which Is reduced by an additional dis- 
placement of, say, 3 mm. This will have little Importance to lift at 
cruise but can be more Important at high angles of attack. 

The raindrops fall on an established film and cause an effective roughness 
as does the existence of waviness In the downstream flow. So far, all 
theoretical attempts have made use of an equivalent sand-grain roughness 
and Table 2, taken from Ref. 4, shows the sand grain roughness equivalent 
for a range of rainfall. The corresponding Increase In skin-friction drag 
Is shown on Table 3 and Is appreciable, although unlikely to be Important In 
terms of fuel consumption for the assumed limited period of the heavy 

rainfall. The variations In maximum lift coefficient and stall angle 

associated with the two forms of roughness are shown In Table 4. Taken 

together with the expected modified airfoil shape, the effects of the 
rainfall are clearly Important at high angles of attack. 

The experimental evidence of the effect of rain on airfoil performance Is 
meager and sometimes contradictory. For example, the tests of Ref. 2 on 
an NACA 64-210 airfoil at R c = 2.6 x 10 & with slat and flaps extended 
showed that a lift loss of up to 30% was possible but results with the same 
configuration at R^ = 1.8 x 10 showed a much smaller lift loss. The 

reason for this difference Is not known but the possibility of the runback 
water clogging the flap gaps had been mentioned as a possibility. It has 
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Table 1. Average film thickness for a symmetric airfoil and 
fuselage at 0 deg angle of attack, 10 m chord 


Rainfall 
rate , 
mm/h 

Calculated 

thickness 

airfoil, 

mm 

estimated 

thickness 

fuselage 

mm 

100 

<0.2 

<0.2 

200 

0.5 or less 

0.2 or less 

500 

0.8 

0.6 

1000 

1 .0 

0.9 

2000 

1 .3 

1.1 


Table 2. 

Equivalent sand-grain roughnesses by rainfall rate 
on a wing 

Rain rate 
mm/h 

k c , mm 

Drop 

Impact 

Cratering 

Waviness 

100 

0.18 

<0.3 

200 

0.37 

0.7 

500 

0.89 

1.2 

1000 

1.83 

1.5 

2000 

3.65 

2.0 


Table 3. 

Increase In total drag due 
friction drag (747 aircraft 

to Increased wing and fuselage 
landing configuration) 


AC 

D /C D 0 ’ X 


Drop 


Rainfall 

Impact 


rate 

Cratering 

Waviness 

100 

1.6 

2.1 

200 

2.3 

3.2 

500 

3.5 

3.8 

1000 

4.6 

4.2 

2000 

5.9 

4.6 
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Table 4. Reduction In maximum lift coefficient and 
angle of attack at stall due to roughness 



SC L 

/c L , % 

Aot r 

t 

max 

deg 

Rain rate. 

Drop Impact 

F 1 lm 

Drop Impact 

FI lm 

mm/h 

cratering 

waviness 

cratering 

waviness 


100 

7 

11 

1-2 

1-3 

200 

13 

20 

1-3 

2-4 

500 

25 

25 

2-5 

2-5 

1000 

29 

28 

3-5 

3-5 

2000 

34 

30 

3-6 

3-5 


also been pointed out that wind tunnel experiments simulating flight in 
rain should be properly scaled In order to model full-scale conditions and 
that this Involves careful consideration of the transition process and of 
wind-tunnel characteristics. The added Influence of the type of surface 
has been demonstrated In Reference 3 for a low Reynolds-number airfoil 
(Wortmann FX-67-K-170) with the resulting lift and drag coefficients of 
Figure 1. In these experiments, the equivalent rainfall was 440 mm/h and 
the chord Reynolds number 310,000. It Is evident from the figure that the 
maximum lift and minimum drag are obtained with a dry surface and that the 
combination of simulated rain and a range of surface coatings Is to reduce 
lift and Increase drag. The surface with a clean epoxy gel may be 
regarded as closest to that of a commercial aircraft but the addition of 
wetting agents Is relevant to surfaces which have been deiced or washed 
with detergent. It Is particularly Important to observe the magnitude of 
the decrease In lift, which occurs close to the angle of attack corres- 
ponding to maximum lift. 

The effects of Figure 1 are. In some measure, particular to the low 
Reynolds number of the Investigation so that It may be that the location 
of transition has been moved forward by the simulated rain. Experiments 
were performed with three different positions of a boundary-layer trip and 
led to the results of Figure 2, which shows that effects of similar mag- 
nitude to those of Figure 1 can be achieved In the absence of the simu- 
lated rain by tripping the boundary layer at locations up to midchord. 
The trip was a 2mm-w1de strip of sand grains of 0.3mm average size so 
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Fig. 1. Effect of surface conditions on lift and drag coefficient as a 
function of angle of attack (Ref. 3). 



Fig. 2. Effect of boundary layer trip on the lift and drag coefficients 
as a function of angle of attack (Ref. 3). 

that, as can be seen from Table 2, It Is similar to the equivalent sand- 
grain roughness of the simulated rain. The position of transition on the 
wings of commercial aircraft Is usually close to the leading edge, so that 
the nature of the surface Is less Important but the differences between 
the dry results and those obtained with simulated rain and the gel coat, 
Figure 1, Indicate the likely effect of heavy rain. In addition, novel 
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designs Involving procedures to maintain laminar flow must take account of 
the Implications of figures 1 and 2. 

Experimentally based Information of the effects of Insect contamination Is 
less than that for rain or Ice, Is confined to low speeds and has been 
considered mainly In relation to the location of transition. This Is an 
Important aspect of laminar control since Insect contamination first 
appears near the leading edge so that laminar flow can be lost if a crit- 
ical level Is exceeded. It has been found that Insect contamination acts 
as distributed roughness and In this sense It may be treated by computa- 
tional methods once the statistical properties of Insects are known. 
Experiments Indicate that, for contamination to occur, a minimum speed, 
probably different for different species, has to be exceeded so that the 
Insects burst on Impact and adhere to the surface. Impact regions or 
capture areas can be calculated by computer programs used for water- 
droplet trajectory and Impingement calculation by substituting for Insect 

mass and drag coefficient. Experiments made on 5-foot-chord airfoils at 
6 

R c ~ 7 x 10 with fruit flies showed maximum roughness heights In 
the range of 0.015 0.030 Inches can occur near the stagnation point 

(Ref. 5). It was found that the roughness height decreased rapidly with 
distance from the leading edge and the effect of this local buildup on the 
airfoil performance was small. The experiments were, however, restricted 
to low angles of attack and larger effects can be expected at large angles 
because the area of peak velocity may occur not far from that of Insect 
Impact. Quantification of this hypothesis can be obtained with the calcu 
latlon method of Section 4. Although no direct tests of Insect contami- 
nation effect on maximum lift are available. It can be assumed that It 
would be similar to that of a distributed roughness. 

It Is evident from Refs. 6 to 10, that Ice accretion may affect the aero- 
dynamic characteristics of airfoils by reducing C. and Increasing 

drag. Two types of Ice may form, rime Ice where low temperatures and 
velocities allow supercooled water droplets to freeze on Impact with a 
resulting accretion similar to that of Figure 3a and glaze Ice at temper- 
atures just below freezing so that water droplets flow along the surface 
and freeze to glaze-ice forms similar to those of Figure 3b. Both types 
of Ice can Influence lift and drag considerably by the modified shape of 
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Fig. 3. Typical (a) rime and (b) glaze Ice accretions on the leading-edge 
of an al rfol 1 . 

the effective airfoil; they also contribute to added drag, and reduced 
lift, through surface roughness. 

The Influence of Ice accretion on the aerodynamic properties of airfoils 
has been Investigated experimentally over many years. In the 1950's, the 
NACA Investigated the effects of Ice on airfoil performance and some of 
these results can be found In the work of Gray and Von Glahn [6] who 
examined NACA 65-212 and 65A004 airfoils and showed the adverse effects of 
Ice on the Integrated lift, drag and moment. Other Investigations, such 
as that by Korkan et al. [7], quantified the effects of simulated Ice 
shapes on airfoil performance and, more recently, Bragg and Colrler [9] 
simulated a measured glaze Ice accretion on a wooden 21-Inch chord NACA 
0012 airfoil and reported measurements of surface pressures, lift and 
moment coefficients and a wake survey to provide airfoil drag. The sepa- 
ration bubble was explored by measuring the tlme-averageo velocities using 
a split-film probe and velocity profiles were obtained In the separation 
bubble for several chordwlse stations at three angles of attack. The 
results show that the Ice shape caused a severe reduction In lift and sub- 
stantial Increase In drag. 

A comprehensive Investigation of the efforts of accretion of frost and 
various Ice formations has been reported by Roed [10] who presents varia- 
tions In lift coefficient measured with a single airfoil, with a trailing 
flap and with extended slat and flap. The measurements [11] were obtained 
In a flight test and show very large modifications of the curve of lift 
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against angle of attack, particularly with the trailing edge extended. The 
angle of attack corresponding to C. can be changed from 9 to 1 degree 

J I In <1 X 

with corresponding reduction In C from 3.5 to 2.8. With the multi 

L me ix 

element configuration, corresponding reductions In C L ^ from 4.3 to 3.3 
were observed. 

It Is evident that the main effect of Ice accretion Is to change the shape 
of the airfoil and so modify Its performance. The prediction of the flow 
characteristics which result from the accretions can be achieved by the 
solution of Invlscld and viscous- flow equations and Interaction of the 
solutions or of the Navler- Stokes equations. As has been shown In various 
papers, for example Reference 12, Interaction between the Invlscld and 
viscous flow equations becomes Increasingly necessary as the angle of 
attack Is Increased. In addition. It Is desirable to make provision for 
the roughness of the Ice surface In a general manner which will also 
accommodate the related roughness effects of rain and Insect contamination. 
To provide this generality, and to permit the Inclusion of accretion model, 
an Interactive boundary layer procedure, based on the solution of differ- 
ential equations In finite- difference form Is advocated and a preferred 
approach Is described In Section 4. 

3.0 Data Correlations 

Roughness caused by rain, frost, snow or freezing fog adhering to the wing 
surface, large accumulations of Insect debris and badly chipped paint can 
play an Important role on aircraft flight performance. These adverse 
effects are addressed In the Federal Air Regulations and have received 
considerable attention In the past several years. Due to the Immense com- 
plexity of the problem, however, estimation of the roughness effects are 
presently limited to data correlations. Computational methods which offer 
broader applicability, accuracy an fundamental understanding are very new, 
as discussed by Shaw [13,14] and their development has so far been limited 
to airfoils. Before we discuss these recent and advanced computational 
methods for airfoils and their possible extension to wings, empennage, 
propellers, rotors and eventually for complete aircraft configurations. It 
Is useful to review the correlations which provide Insight Into the effects 
of small amounts of wing-surface roughness on aircraft flight performance. 
In addition, the shortcomings of correlations for predicting the effects 
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ol Ice ori lift and drag of airfoils are considered. Previous reviews of 
performance degradation of propellers, helicopter rotors (hover and forward 
flight) and complete aircraft are available In References 14 to 20. 

As discussed by Brumby [21] for full wing-span upper surface roughness 
beginning at the leading edge and extending varying distances aft, typical 
effects are a reduction of the maximum lift coefficient (Increase In stall 
speed), a reduction In the angle of attack at which stall occurs and a 
rapid post-stall drag Increase (see Fig. 4). The effects become more 
severe as the size and chordwlse extent of the roughness Increase and they 
may be accompanied by a reduction In lift at a given angle of attack and 
by an Increase In the parasite drag. 

Figure 5 shows Brumby's correlation of wind tunnel and flight data and the 
effects of surface roughness on the maximum lift coefficient of a wing. 
The majority of the data are from two-dimensional tests but the four 
flagged points represent data obtained from three-dimensional swept sur- 
faces and appear to confirm that the correlation Is applicable to wings. 
The data are for two general types of roughness on wings without leading- 
edge high- lift devices. The solid symbols Indicate data for distributed 
(sand-grain type) roughness at the leading edge, or on the entire upper 
surface, and the open symbols correspond to localized full-span 




Fig. 4. Typical effect of surface roughness at the leading edgeon aero- 
dynamic characteristics. 
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Fig. 5. Reduction of maximum lift coefficient due to wing surface 
roughness . 


TYPE OF ROUGHNESS 


SAND GRAIN BAND 
PROTRUDING STRIP 
MULTIPLE GROOVES 
CARBORUNDUM GRIT 
SAND GRAINS 
WIRE MESH ON SURFACE 
FWD FACING STEP 
PROTRUOING STRIP 
PROTRUDING STRIP 
PROTRUDING STRIP 
CARBORUNDUM GRIT 
FROST (IN ICING TUNNEL) 
INSECT CONTAMINATION 
SIMULATED TAIL PLANE ICE 
simulated TAILPLANE ICE 
simulated deicer boot 

CARBORUNDUM GRIT 
CHIPPED PAINT ON IE 
BALLOTINI 

BURRED RIVETS ON LE 

simulated frost 

SIMULATED WING ICE 

simulated ice roughness 


disturbances at various chordwlse stations. It Is clear that large 

decreases In C, , with resultant large Increases In stall speed, can 
L max 

occur due to comparatively minor wing surface disturbances. 


Bragg et al. [20] claim that Brumby's correlation, while useful In estima- 
tion of changes In C, , Is limited In that It contains no Reynolds 
number effects and little detail of the actual roughness or Its density. 
As a consequence, It falls to predict the measured results of Increased 

C, due to Ice accretions. 

L max 

An early correlation equation for the effects of Ice on drag was formulated 
by Gray [22] and based on the data collected In the NASA Lewis Icing 
Research Tunnel (1RT) primarily In the 1950's. It Is best suited for glaze 
conditions and Fig. 6, taken from Ref. 14, shows that It provides guidance, 
though the predicted drag rise Is often too large. More recently, Bragg 
[23] developed a rime-ice correlation, also based on the data gathered In 
the NACA IRT. Flemming [18] acquired a large data base In the Canadian 
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Fig. 6. Gray's drag rise correlation for two airfoils: (a) NACA 0012 (b) 
NASA MS{ 1 ) -031 7 , both taken from Ref. 14. 


NRC high-speed Icing wind tunnel for a series of reduced scale rotor air- 
foil sections which he used to develop a series of correlations for the 
drag rise due to Ice accretions. Figure 7 shows Flemming's correlation for 
the same two airfoils as Fig. 6, also taken from Ref. 14 and that contrary 
to Gray's correlation, It provides a lower drag rise. 

It Is evident from the results discussed above and from the additional 
examples of Ref. 14, that the experimental data Involve effects not repre- 
sented by the correlation equations. They are confined mainly to drag and 
do not Include terms to take account of known effects such as those of 
Reynolds number, airfoil shape and slats. There Is a clear need for a 
procedure which will represent the aerodynamic properties of the flow 
around airfoils correctly and will allow correct representation of large 
changes In geometry, such as those associated with accretions of rime and 
glaze Ice, as well as the smaller changes associated with frost, Insects 
and rain. It Is also desirable that this procedure should be able to deal 
with the three-dimensional effects of real airplanes. The following sec- 
tion addresses these needs. 
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Fig. 7. Flemming's drag rise correlation for two airfoils: (a) NACA 0012 
(b) NASA MS( 1 ) - 031 7 , both taken from Ref. 14. 

4.0 Numerical Solution Procedures 

The development of digital computers, and their ability to solve many alge- 
braic equations In a short time, has spawned considerable efforts to solve 
the conservation equations of fluid mechanics and heat transfer In differ- 
ential form. In the field of aerodynamics, these numerical solution 
procedures have been directed to the solution of reduced forms of the 
Navler Stokes equations and particularly to Interactive boundary- layer and 
thin layer Navler-Stokes equations. Reviews, such as that of Cebecl and 
Whltelaw [25J show that useful calculations can already be performed for 
two- and three-dimensional flows and over an extensive range of angle of 
attack . 

Two approaches have been used to obtain the results of Section 5, and are 
based on the Interactive boundary layer procedure of Cebecl et al. [12] and 
on the thin layer Navler-Stokes procedure developed at NASA Ames [26, 27]. 
The following paragraphs provide outlines of these two approaches and Sub- 
section 4.1 a more detailed description of the Interactive procedure, which 
has been used to obtain most of the results of Section 5. The turbulence 
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model and the modifications necessary to deal with rough surfaces are 
considered In Subsection 4 .2 . 

Ihe Interactive procedure Involves solution of Invlscld flow equations and 
of boundary layer equations. The results of Section 5 were obtained with 
steady, two dimensional equations. The Invlscld flow was determined by 
conformal mapping and by a panel method and the boundary-layer equations, 
with turbulent diffusion represented by an eddy viscosity approach, were 
solved by the two-point finite-difference method of Keller. Interaction 
between the Invlscld and viscous flows was achieved by a blowing-velocity 
distribution which was linked to the displacement- thickness distribution 
through the Hilbert Integral. Where separation was encountered, the equa- 
tions were solved In Inverse form with the FLARF. approximation which neg- 
lects longitudinal convection In the recirculation region. The same 
approach was taken In the wake where the dividing streamline was computed 
from the Invlscld flow as a line having constant stream function and was 
used by the Invlscld method to apply the blowing velocity required to 
simulate the displacement thickness and to compute the Invlscld velocity. 
An Invlscld point distribution In the wake was defined on which the wake 
velocity distribution was determined. This required Interpolation of the 
boundary- layer blowing velocity onto the Invlscld points. In addition, the 
computed Invlscld velocity was Interpolated back to the boundary- layer 
points by making use of the computed velocity at the trailing edge as the 
Initial wake point. In the Immediate vicinity of the trailing edge, par- 
ticular care was required In the choice of the locations at which values 
of the blowing velocities were applied In the solution of the Invlscld 
equations. Further details are provided In the following subsection. 

The thin layer Navler Stokes equations are generally referred to In the 
literature as TLNS equations and are obtained by neglecting streamwlse and 
spanwlse derivatives of the viscous and turbulence stress, conductive heat- 
flux terms, and any term Involving mixed derivatives. These approximations 
are Justified either by order of magnitude arguments or by consideration 
of computational accuracy argument [28,29]. The TLNS equations have been 
proposed mainly on the computational argument. The form of the equations 
generally used does not satisfy relationships between metric coefficients 
In diffusive and conduction terms but the resulting error Is usually 
Insignificant, except when the effective viscosity Is relatively large. 
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In addition, longitudinal- curvature diffusive terms are neglected as a 
consequence of the Cartesian velocity components. 

An Implicit numerical method was used to solve the 1LNS equations and Is 
based on that of Ref. 30. Since only steady-state computations are of 
Interest, a diagonal form for the Euler equations and a spatially varying 
time step were used. Reference 27 provides a description of the numerical 
scheme to obtain the results presented here. The turbulence model Incor- 
porated In the code Is described In Ref. 31. 

4 . 1 Interactive Boundary-Layer Procedure 

The outline of the Interactive approach provided above Is expanded In this 
section so as to allow more detailed assessment of Its features. The 
conformal-mapping method for the solution of the Invlscld-f low equations 
has been described extensively In the paper by Halsey [32] and the panel 
method by Hess and Smith [33] to which the reader Is referred for further 
Information. 


The boundary layer equations are solved In their two-dimensional form, 

3u dv _ 

T7 f ^ - 0 

3x 3y 


( 1 ) 


u 


3u 

3x 



du 


e dx 


3 3u. 

3y <b P 


( 2 ) 


where b = 1 * e /v and c Is defined by a form of the eddy viscosity formu- 
m m 

latlon of Cebecl and Smith [34] discussed In Subsection 4.2. For wall 
boundary layer flows, the boundary condition may be written as 


y - 0 , u = 0, v = 0; y > 4, u > u e (3) 

and for the asymmetric wakes of airfoils, 

y -> - 64 , u = u e , y = y*. v = 0, y » 4 U , u = u e (4) 

Here 4 and 4 denote the lower and upper wake boundary- layer thicknesses, 
8 . u 

respectively, with y* representing the dividing streamline assumed to be 
given. The above equations assume that there Is no pressure gradient 
across the shear layer but the corresponding constraint can be removed for 
strongly curved wakes. 
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The solution of Eqs. (I) and (2) can also be obtained by a procedure In 
which the external velocity Is computed as part of the solution. This 
procedure Is known as the I nverse problem and Is essential to remove the 
singularity associated with an external boundary condition based on a 
specified distribution of u g . It Is necessary to specify an additional 
boundary condition In addition to the boundary conditions given by Eqs. (3) 
and (4) since u (x) now represents an unknown and. In the Interactive 
scheme of Cebecl et al. [12], this Is accomplished by rewriting the exter- 
nal velocity, u g (x) as 


u e (x) = ug(x) + 6u e (x) 


(5) 


where u°(x) denotes the Invlscld velocity and 6u(x) Is the perturbation 

g' ' c 

velocity due to viscous effects. The latter Is related to the blowing 
velocity Induced by the boundary layer by a variation of the Hilbert 
Integral 


Su g (x) 




**) 


do 

s o 


( 6 ) 


with the Interaction region limited to a finite range x fl < x < x fa . 


Following Cebecl and Clark [35], we write Eqs. (5) and (6) as 


U e (x) . U°(x) * c 1J (u e 4*) J 


(7a) 


Here c denotes the Interaction-coefficient matrix, which Is obtained 
from a discrete approximation to the Hilbert Integral In Eq. (6). In this 
form, Eq. (7a) represents an outer boundary condition for the Inverse prob- 
lem. It can be generalized to the form 


U e (x) • u*(x) * ^ - (».•*),!. 


(7b) 


where u (x) corresponds to the Invlscld velocity distribution which 
® k 

contains the displacement effect (6*) computed from a previous sweep. 


The solution of Eqs. (1) and (2) for the Inverse problem Is now obtained 
subject to Eqs. (3) for wall boundary layers and to Eqs. (4) for wake flows 
with u g (x) given by Eq. (5). It Is convenient to express the above 
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equations In terms of transformed coordinates and, with u denotlnq a 

o 3 

constant reference velocity, we Introduce the transformation 

u 1/2 

n = y. ♦ = (u Q vx) /2 f(x,n) (8) 

With primes denoting differentiation with respect to n, Eqs. (1) and (2) 
and their boundary conditions on the airfoil and In the wake can be written 
In the following form: 




On the airfoil 
n = 0, 


n = V 


In the wake 
i = - n a 


n 

where 


V 


f = f 1 =0 
f' = w, 

f' = w; 
f' - w. 


w - c u (n e v - f) = 9l 

n = n*, f = 0 

w " c 11 [w(n u ~ V “ < f u " = 9 1 


’ll 


c ( ^) 1/2 
C 11 u ] 


( 8 ) 

(9a) 

(9b) 

(10a) 

(10b) 


Here w denotes the dimensionless external velocity u /u and the par- 

e o 

ameter g^, which results from the discrete approximation to the Hilbert 
Integral Eq. (6), Is given by 


k 1 - 1 k 

9i - w ♦ l c ^ j ( D j - DJ) - c^D 


J = 1 


•iri 


where 


n /Mv 1/2 , , . 

D = <u > ( V - V 
0 


( 11 ) 


( 12 ) 


The expression for g^ on the wake Is nearly identical to that for the 
airfoil, Eq. (11), except that now Eq. (12) Is given by 


D = C> 1 / 2 [w(n u V (f u - f % )) 


(13) 


The numerical solution of Eq. (8) subject to the boundary conditions given 
by Eqs. (9) and (10) has been obtained with Keller's Box method which Is 
an efficient, second-order finite-difference method extensively used by 
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Cebecl and his associates for a wide range of flows, as discussed In [36]. 
The procedure for wake flows Is novel and has some desirable features which 
allow the calculation of thick separated boundary layers by the extension 
of the Interactive scheme. As In the solution of wall boundary-layer 
flows, we assume 


t = u (14a) 

u ' = v (14b) 

and write Eq. (8) as a first-order system 

.... 1 , dw „ 3u 3f 

(bv) * ~ fv r XW - d - , X(u ^ v ^) (14c) 

The Mechul function formulation [37] Is used to obtain stable solutions to 

the above equations and to reduce their sensitivity to the boundary condi- 

tions which Involve and f . Since both s and w are functions of 
w only, we write 


s ’ = 0 ( 14d) 

w' = 0 (14e) 

and the boundary conditions for the system given by Eq. (14) become 

n = -n a , u = w, s = f s. ; n = n*. f = 0 (15a) 

n = V u = w, w - c^ [w(n u - n^) - (f - s) ] = g 1 (15b) 


The system of Eqs. (14) and (15) are solved by the procedure described In 
Ref. 36. After the finite- difference approximations to Eqs. (14) are 
written, the resulting nonlinear algebraic system Is linearized by Newton's 
method and the linear system Is then solved by the block-elimination 
method. 

In computing airfoil flows with separation [12], again the FLARE approxi- 
mation was used to obtain stable solutions on the airfoil and In the wake. 
As the extent of the separation region Increased, however, Cebecl et al. 
[12] Introduced an additional Iterative scheme based on a continuation 
method at the start of the wake calculations. With u ^ corresponding 
to a nonseparating velocity profile constructed somewhat arbitrarily from 
the separated velocity profile at the trailing edge an Initial velocity 
profile was defined by 
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u = u ref <■ n(u u ref ) n -= 0, 0.50, 1.0 (16) 

and the boundary- layer solution was computed at the first point on the wake 
with n = 0. The solution was repeated for other values of n until converg- 
ence. The procedure was applied for each profile In the wake with separa- 
tion and was necessary at higher angles of attack near stall conditions. 
For additional details, see Choi [38]. 


4-2 Eddy-Viscosity Formulation 

The presence of In b requires a turbulence model and In Ref. 12, the 
algebraic eddy- viscosity formulation of Cebecl and Smith [34] was used for 
"clean" airfoils. According to this formulation for wall boundary- layer 
flows, Is defined by two separate formulas, given by 


ft 2 


dU 

ay 


'tr 




a 


J (u o - u)dy 


T tr T 


o < y < y r 


y c - y - 4 


where with k = 0.4 


(17a) 


(17b) 


L = icy [1 - exp( - y/A) ] 


A = 26uu 


-1 


u = (— ) 1/2 f 
t p 'max 


3u 

= p 3y ' 


1 


1 + 5. 5(y/6) 


(18) 


The condition used to define y^ Is the continuity of the eddy viscosity; 
from the wall outward Eq. (17a) Is applied until Its value Is equal to the 
one given by Eq. (17b). In Eq. (17), Y^ r Is an Intermlttency factor 
which accounts for the transitional region between a laminar and turbulent 
flow and Is given by 


Y tr = 1 - exp[-G( x 


V 


I ^] 

u J 
x tr e 


(19a) 


Here G Is an empirical parameter which, with x denoting the location 
of the start of transition and R x^ r the transition Reynolds number 
defined by R x^ r = (u^/v)^, Is given by 
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J 

1200 


R 


- 1.34 
( tr 


(19b) 


According to the Cebecl-Smlth (CS) model, the parameter a In Eq. (17b) 

Is equal to 0.0168 for values of R. greater than 6000, and Is given by 

the expression In [34] for R less than 5000. Studies Indicate, how- 

0 

ever, that In flows with strong pressure gradient, the value of a should 
also be changed when R Q > 5000. For this purpose Cebecl et al. use an 
expression for a to account for strong adverse pressure gradient effects 
as discussed In Ref. 12. 


The above eddy-viscosity formulation of Cebecl and Smith for "clean" air- 
foils can also be used for "rough" airfoils with small modifications to the 
Inner eddy- viscosity formula of Eq. (17a) only. We use the formulation of 
Cebecl and Chang [38] for this purpose and rewrite the modified mixing 
length expression of Eq. (18) as 


L = *(y *■ Ay ) { 1 exp[-(y «■ Ay)/A] } 


( 20 ) 


where Ay Is a function of an equivalent sand-grain roughness k . In 

terms of dimensionless quantities with k* = k u /v, we have 

s s T 


f0.9 [ - k + exp (- k + /6) ] 5 < k* < 70 


Ayu 


f 0.58 
0.7 (k $ ) 


70 < k < 2000 
- s - 


(21) 


5.0 Results and Discussion 


The calculated results are presented In three subsections which deal with 
smooth airfoils, rough airfoils and Iced airfoils, respectively. The first 
subsection Is Included to quantify the extent to which the two calculation 
methods can represent airfoil flows as a function of angle of attack and 
without the added complication of a roughened surface. The rough surfaces 
of section two allow examination of the value of the concept of equivalent 
sand- grain roughness within the framework of the turbulence model of the 
Interactive procedure. Similar results can be expected from solution of 
the thin-layer equations and from the application of both procedures to 
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problems of Iced airfoils were the Ice can be considered as roughness. The 
last section describes results obtained with the two calculation methods 
for Iced airfoils where the Ice accretion changes the shape of the leading 
edge of the alrfol 1 . 

5 . 1 Smooth Airfoils 

Figures 8, 9 and 10 present results obtained with the Interactive and thin- 
layer procedures for two airfoils and angles of attack up to around 16 
degrees. They are taken from Ref. 39 In which additional results can be 
obtained . 


Measurements of the flow around a NACA 4412 airfoil have been reported In 
Refs. 40 and 41 and made use of flying hot wire anemometry at angles of 
attack up to that of maximum lift. The Reynolds number based on chord 



Fig. 8. Variation of lift coef- 
ficient with angle of 
attack for the NACA 4412 
airfoil at R c .•= 1 .523 
x 1 06 . 




Fig. 9. Variation of total drag coef- 
ficient with angle of attack, 
(a) R c = 1.523 x 10 6 , (b) 

R c = 3 x 106. 
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Fig. 10. Variation of lift and drag coefficients for these GA(W)-2 air- 
foil, R c = 4.3 x 10&. 

6 

length was 1.523 x 10 and transition was Induced at 2.5% chord on the 
upper surface and 10.3% chord on the lower surface. The same Reynolds num- 
ber and transition locations were used In the calculations with both meth- 
ods at angles of attack up to 12 degrees. Above 12 degrees, the Inter- 
active calculations revealed a laminar separation very close to the leading 
edge and the location of the onset of this separation was taken as that of 
the onset of transition. This assumption Is consistent with transition 
having occurred upstream of the trip In the experiment. 


The flow around a NACA 4412 airfoil has also been Investigated at a chord 
Reynolds number of 3 x 10^ (Ref. 42) and corresponding calculations have 
been performed with the onset of transition determined. In the absence of 
experimental Information, by Michel's formula [43] given In Ref. 44, that 
Is 


R 0 = 1.174 (1 


22.000 . d 0.46 

R ' x 
x 


( 22 ) 


or by the onset of laminar separation. 


Figures 8 and 9 permit comparison of the measured and calculated values of 
lift and drag coefficients for the NACA 4412 airfoil. The results for lift 
coefficient display the variation with angle of attack obtained from the 
solution of the Invl scld-f low equations and which diverges from the mea- 
surements with Increasing angle. It Is clear that the two calculation 
methods agree well with each other and with experiment up to around 8 
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degrees beyond which the Interaction procedure follows the experimental 
results more closely and represents the expected maximum value at the same 
angle of attack on the measurements. This aspect of the comparison of the 
two calculation methods Is similar to that reported In Ref. 31 In relation 
to a NACA 0012 airfoil, with the two calculation methods agreeing well with 
each other up to around 14 degrees. 

The drag coefficients of Figure 9 are more difficult to appraise since the 
two experimental distributions differ Increasingly with angle of attack 
and, for example, by a factor of almost two at 12 degrees. Experimental 
differences are common, due to the Inaccuracy of the Integrating the wake 
profile and to wind-tunnel effects associated with blockage or finite span. 
The magnitude of the present differences Is, however, unusual. The results 
obtained with the two calculation methods agree well with the measurements 
of Refs. 42 and 45 and are at odds with those of Refs. 40 and 41. They 
were obtained by wake Integration of the Interactive boundary-layer results 
and by surface Integration of the TINS results since. In the latter case, 
the far wake was not well represented by the calculations. 

It Is Interesting to note that measurements of drag coefficient obtained 
with NACA 0012 airfoils and described In Refs. 46 and 47 also show d 1 s 
crepancles which Increase with angle of attack so that, at 12 degrees, the 
difference was around 15%. The calculations of Ref. 31 revealed similar 

discrepancies and It Is clear that It Is difficult to achieve a high degree 
of accuracy. 

Distributions of pressure coefficients agreed very closely at 0 and 4 
degrees angle of attack and agree reasonably well at 12 degrees but with 
differences In the trailing edge region and particularly In the repre- 
sentation of the small region of upper surface separation which affected 
the wake and led. In part, to better representation of the wake by the 
Interaction boundary-layer method. It was found that the pressure peak was 
very sharp and located almost at the leading edge. The distributions at 
16 degrees revealed an even more peaky pressure distribution and the need 
to specify the onset of transition upstream of the trip can be appreciated. 
The discrepancy between the two calculation methods was more evident at 16 
degrees and Involved a large and Important difference In the tral 1 Ing-edge 
region where the Interactive approach suggested separation some 20% of 
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chord upstream of the trailing edge and the TLNS method suggested a smaller 
separation. 

The GA(W) 2 airfoil represents a more difficult test since It Is a 13% 
thick airfoil of supercritical form. Measurements have been reported [48] 
at a chord Reynolds number of 4.3 x If)* 5 and with transition trips at 
7.5% chord on both surfaces. Calculations were performed with both methods 
and, as before, considered the onset of transition In accord with the 
experimental configuration unless the Interactive approach Indicated a 
laminar separation bubble. In which case the onset of transition was taken 
as coincident with the onset of separation. The results shown In Figure 10 
display the variations of lift and drag coefficient with angle of attack. 

The deductions which can be made from Figure 10 are similar to those 
obtained from Figures 8 and 9. The two calculation methods represent an 
Improvement over the Invlscld-f low calculations In terms of the lift coef- 
ficient and, as before, the Interactive method provides results In very 
close agreement with experiment. The Interactive method Is also able to 
calculate the drag coefficient with accuracy which diminishes with angle 
of attack so that the difference between measurement and calculation Is 
around 20% at 14 degrees. The results of the TLNS method are less satis- 
factory. A sample of pressure distributions Is provided on Figure 11 and, 
as before, shows the need for proper Implementation of the TLNS method In 
the tral 1 Ing- edge region. 

Before we conclude the discussion on clean airfoils. It Is useful to point 
out the Importance of Including the wake effect In the Interactive 
boundary- layer method. Studies by Cebecl et al. [12] Indicate that In 
high Reynolds number flows over airfoils at small and moderate angles of 
attack. It Is sufficient to perform the calculations up to the trailing 
edge. The effect of the wake becomes Important at higher angles of attack, 
especially In flow conditions approaching stall angle, and must be 
accounted for In the calculations. As can be seen from the results shown 
In Figure 12a for the NACA 0012 airfoil, the effect of wake on the 
tral 1 Ing-edge displacement thickness Is negligible for a = 10° but more 
than 30% for a = 16° Indicating that without the wake effect, the mag- 
nitude of the tral 1 Ing-edge displacement thickness Is significantly greater 
than Its value with the wake effect. The reduction of the displacement 
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Fig. 11. Comparison of calculated and experimental pressure distributions 
for the GA(W) -2 airfoil, R c = 4.3 x 10 6 . (a) o = 0°, (b) 

a = 6°, (c) a = 12°. 

thickness reduces the flow separation on the airfoil and decreases the lift 
coefficient, as shown In Figures 12b and 12c, respectively. 

The effect of wake Is also Important at low Reynolds number flows, which 
are dominated by large regions of separation bubbles leading to relatively 
large tral 1 Ing-edge displacement thicknesses even at low angles of attack. 
This situation Is analogous to high Reynolds number flows over airfoils at 
high angles of attack and again It requires the Inclusion of the wake In 
the calculations. Further details are provided In Ref. 49. 

5 . 2 Rough Airfoils 

To examine the ability to deal with airfoils with rough surfaces, the 
experiments of [50] were represented by the Interactive boundary-layer 
method. The roughness comprised 0.001 Inch carborundum grains applied to 
24-Inch chord airfoils and spread evenly over a surface length of 0.08 
chord. Within the framework of the turbulence model of Subsection 4.2, It 
Is necessary to convert this form of roughness Into equivalent sand-grain 
roughness. This was done with the procedure of Smith and Kaups [51] In 
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Fig. 12. Effect of wake on the (a) displacement thickness, (b) separation 
region, and (c) lift coefficient for the NACA 0012 airfoil, R c 
= 6 x 106. 


which the ratio of the equivalent sand grain roughness to the roughness of 
the applied elements, k $ /k, was assumed to be a function of the concen- 
tration and shape of the roughness elements, see Fig. 13. In all cases 
considered here, the shape of the elements was approximated by a sphere, 
and the concentration, which represents the mean value of the area covered 
by the roughness elements was taken as 0.076 and the equivalent sand-grain 
roughness height was obtained from Fig. 13 to be 

k s /c = 0.00094 (23) 

The above expression was used for three airfoils NACA 4412, 23012 and 0012 
for which experiments were performed at a Reynolds number of R = 6 x 
10 6 . 
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Figures 14 and 15 contain the results 
of the computations for the NACA 4412 
airfoil. It can be seen from the 
lift coefficient results of Fig. 14a 
that the computations and measure- 
ments agree well up to the stall 
angles for the smooth airfoil as well 
as for the airfoil with leading-edge 
roughness. The drag curves of Fig. 
14b also show good agreement between 
computations and measurements 
although for the clean airfoil the 
drag Is slightly underpredicted at 
higher angles of attack. Values of 
dimensionless displacement thickness 
A*/c at the trailing edge are 
shown In Fig. 15a and, since the 
transition location Is much further 



Fig. 13. Equivalent sand-grain rough- 
ness for uniform three- 
dimensional roughness as a 
function of concentration. 
Dashed lines are extrapola- 
tions of experimental data 
[51]. 




Fig. 14. Comparison of calculated and experimental results for the NACA 
4412 airfoil at R c = 6 x 10&. 
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Fig. 15. Results for the NACA 4412 airfoil at R c , 6 x 10 6 (a) trail- 
ing edge displacement thickness, (b) extent of flow separation on 
functions of angle of attack. 


upstream In the case of a rough airfoil than In the case of a smooth air- 
foil, the trailing edge value Is higher In the rough case with C 
consequently lowered. Figure 15b shows the extent of the tral 1 Ing-edge 
separation as a function of the angle of attack and quantifies the earlier 
separation associated with the rough surface. 

The NACA 23012 airfoil, like the NACA 4412, has camber and the results of 
the calculations are shown In Figs. 16 and 17. The lift and drag curves 
of Fig. 16 again demonstrate good agreement between the computations and 
measurements, and the drag curve Is In even better agreement than that for 
the NACA 4412. Figure 17a shows the value of dimensionless displacement 
thickness at the trailing edge and Fig. T 7 b the trailing edge separation 
both as a function of the angle of attack. The slope of the tral 1 Ing-edge 
separation curve of Fig. 17b at 11° angle of attack suggests that, with the 
rough surface, separation would occur over a very large portion of the 

airfoil If the angle of attack were Increased. 
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Fig* 16. Comparison of calculated and experimental results for the NACA 
23012 airfoil at R c = 6 x 10&. 
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Fig. 17. Results for the NACA 23012 airfoil at R c = 6 x 10 6 . (a) 

Trailing edge displacement, and (b) extent of flow separation as 
functions of angle of attack. 



Ihese results, together with the results shown In tig. 18 for the NACA 0012 
airfoil, quantify the extent to which the Interactive boundary-layer method 
can represent the flows over airfoils with roughened surfaces. It Is evi- 
dent that the calculation method correctly represents this effect and the 
trends of the lift and drag curves. Where the effects of Ice, rain or 
Insect deposition can be regarded as roughness with an equivalent sand- 
grain value, similar results can be expected. 



(a) (b) 

Fig. 18. Comparison of calculated results for the NACA 0012 airfoil at 
R c = 6 x 10 6 . 

5.3 Iced Airfoil 

This section presents results for airfoils with Ice accretions large enough 
to change the shape of the leading edge. The results of the two calcula- 
tion procedures are compared with the measurements of Bragg and Colrler 
[9], obtained with an NACA 0012 airfoil at a Reynolds number of 1.4 x 10 & 
and at angles of attack up to 10 degrees. A large change In the leading 
edge was arranged with a wooden attachment to represent the shape of a 
typical glaze-ice formation. The calculations with the TLNS equations were 
performed by Potapcynk [52] with the ARC20 code and the grid of 253 x 64 
nodes shown In Fig. 19. 

The lift and drag coefficients computed with the ARC2D code are shown In 
Fig. 20 as functions of angle of attack. As can be seen, the results agree 
well with the measurements up to an angle of attack slightly smaller than 
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Fig. 19- Grid used In the TLNS code 
for the NACA 0012 Iced alr- 
f ol 1 at a = 1 0° . 


that of maximum lift. At this angle, 
as In the case of "clean" airfoils, 
the computed lift coefficient does 
not agree as well but the drag coef- 
ficient agrees remarkably well with 
the data. 

Before we present the results ob- 
tained with the Interactive boundary- 
layer approach. It Is useful to 
comment on the Invlscld method and 
to discuss the role of the Ice on 
the boundary- layer calculations. In 
the latter case the Ice accretion 
can drastically change the pressure 



Fig. 20. Computed (a) lift and (b) drag coefficients for the Iced airfoil 
of Ref. 9. 

distribution near the leading edge and can cause the viscous flow calcula- 
tions to break down. 

The Interactive boundary layer results presented In the previous two sec- 
tions were obtained with Invlscld flow computed by the conformal-mapping 
technique developed by Halsey [32], While this technique gives excellent 
results for clean and rough airfoils, preliminary studies showed that It 
was less satisfactory for Iced airfoils of Interest. The reason appears 
to be that the conformal mapping uses up to 250 points equally spaced 
around the circle Into which the airfoil Is mapped and these are not 
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sufficiently concentrated at the leading edge to represent a complicated 
shape such as that of the accumulated Ice. for this reason, attention was 
directed to a panel method [33] In which the airfoil Is defined by a set 
of points In the physical plane, which would allow the concentration of 
points In the leading edge region and to verify graphically that the Ice 
shape has been adequately represented. Neighboring points on the airfoil 
are connected by straight-line panels so that. In a sense, the airfoil Is 
approximated by a high-order Inscribed polygon, tach panel has both source 
density and vortlclty distributed along It with panel vortlclty strengths 
set equal so that the vortlclty Is defined by a single parameter, total 
strength, which Is adjusted to satisfy the Kutta condition. The source 
strengths, however, have Independent values on each panel and these are 
adjusted, by solving a set of simultaneous linear equations, to satisfy the 
normal- ve loci ty boundary condition at the midpoints of the panels. In the 
strictly Invlscld case this condition requires that the total normal vel- 
ocity, freestream plus body sources and vortices, should vanish. When the 
boundary layer Is simulated, the desired normal velocity Is finite and 
equals the derivative along the surface of the product of tangential vel- 
ocity and displacement thickness. It Is known that this surface blowing 
distribution displaces the dividing streamline outward from the surface of 
the airfoil to the location of the displacement surface. Experience has 
shown that best results are obtained when the surface pressures are calcu- 
lated and the Kutta condition Is applied on the displacement surface, 
rather than on the surface panels. 

In general, boundary layer calculations are rather sensitive to rapid var- 
iations In the external velocity distribution. In order to maintain comp- 
utational accuracy and avoid early breakdown of the solutions In regions 
of steep adverse pressure gradients. It Is necessary to take fine steps In 
the streamwlse direction. For airfoils with large Ice accretions, however, 
It Is further necessary to reduce the sensitivity of the boundary-layer 
calculations to the pressure distribution. In the extension of the Inter- 
active boundary- layer approach of Cebecl et al. to Iced airfoils, this Is 
accomplished by using a continuation method In which the prescribed Ice 
shape Is Introduced Into the calculations gradually. Figure 21 shows a 
sketch of the Iced airfoil In which the Ice shape changes In Increments of 
n ranging from 0 to 1 with n = 0 corresponding to the clean airfoil and 
n = 1 to the airfoil with the prescribed Ice shape of Ref. 9. 
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In performing the Interactive bound 
ary- layer calculations for Iced air 
foils, unlike the procedure used for 
clean and rough airfoils, the viscous 
flow calculations were performed only 
up to the trailing edge and did not 
Include the Influence of the wake, 
thus restricting the accuracy of the 
solutions at high angles of attack. 

Studies are under way to remove this 

restriction. Fig. 21 . Ice shapes used In the con- 

tinuation method of the 
, , , , Interactive boundary-layer 

At first the calculations were per scheme; n = 1 corresponds to 

formed for the clean airfoil (n = 0) the prescribed shape. 

at a = 0. After convergence, the 

Ice shape was Introduced Into the calculations by taking the value of n 
equal to 0.4 and Iterating the solutions until convergence. Subsequent 
calculations were then made for new values of n equal to 0.5, 0.6, 0.7, 
0.8, 0.85, 0.90, 0.925, 0.950, 0.975 and 1.0. Once a complete converged 
solution for a - 0 was obtained, the calculations for another angle of 
attack were performed for n = 1 by Initially computing the pressure dis- 
tribution for the new a for the blowing velocity of the Iced airfoil at 
a = 0. With each solution of the boundary- layer equations, a new blow 
Ing velocity was computed to obtain a new pressure distribution, and, as 
before, this procedure was continued until convergence. At small angles 
of attack. It was sufficient to choose the angle of attack Increments, 
Aa, to be around 0.50°; at higher angles of attack, especially at con- 
ditions approaching stall, Aa had to be chosen smaller, becoming around 
0.1° for a's between 5° and 6°. 


Figure 22 shows the Invlscld external velocity distribution near the lead- 
ing edge of the airfoil shown In Fig. 21. As can be seen, the Invlscld 
velocity distribution differs significantly with and without Ice: the clean 
airfoil has a favorable pressure gradient followed by an almost zero pres- 
sure gradient whereas the Iced airfoil has a severe adverse pressure gra- 
dient after a short Initial region of favorable pressure gradient. For 
both surfaces the rapid flow deceleration Is followed by a gentle favorable 


33 



Mg. 22. Invlscld u e /Uo, distributions for Iced and clean NACA 0012 
alrfol 1 at a = -0. 15°. 

pressure gradient and the locations of possible flow separation can be 
easily Identified. 

The results shown In Figs. 23 to 27 correspond to Interactive boundary- 
layer calculations obtained for the Invlscld velocity distribution of Fig. 
22 after several boundary- layer sweeps along the airfoil. As In the case 
of a clean airfoil, the calculations were started at the stagnation point 
of the airfoil. The leading-edge results of Fig. 23 show that the computed 
external velocity distribution changes drastically with each sweep from 
that predicted by Invlscld flow theory, but the location of flow separation 
for each surface remains essentially unchanged. The calculations, see also 
Fig. 24, Indicate a 10-percent separation bubble for the upper surface and 
a 25-percent separation bubble for the lower surface. Unlike the separa- 
tion point, the reattachment point moves upstream with each sweep, but the 
difference In reattachment points becomes smaller with Increase In the 
number of sweeps. The velocity profiles on the upper surface of the air- 
foil, Fig. 25, show that the extent of flow separation Is large and that 
the present method Is still able to cope well with It. 

Perhaps the biggest surprise In the Interactive flow calculations Is the 
behavior of the displacement thickness distribution on the airfoil. Since 
the Incidence angle Is practically zero and the airfoil Is symmetrical, 
the displacement thickness distributions on both surfaces are the same for 
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Fig. 23. Variation of external velocity distribution on Iced NACA 0012 
airfoil with each boundary- layer sweep for a = -0.15° and R c 
= 1.5 x 10&, (a) upper surface, (b) lower surface. 



Fig. 24. Variation of local skin friction coefficient distribution on Iced 
NACA 0012 airfoil with each boundary-layer sweep for a = -0.15° 
and R c = 1.5 x 1 0 6 , (a) upper surface, (b) lower surface. 

the clean airfoil. Its magnitude at the trailing edge Is approximately 
one-half percent of the airfoil chord, which Is relatively small and has a 
very small effect on the overall pressure distribution. In the case of the 
Iced airfoil, the flow separation due to Ice alters the displacement thick- 
ness distribution on the airfoil, as shown In Fig. 26. The lower surface 
has a large separation bubble which causes the magnitude of the displace- 
ment thickness at the trailing edge to be about two percent of the airfoil 
chord. The upper surface has a smaller separation bubble and, as a result, 
the displacement thickness at the trailing edge Is about one percent of the 
airfoil chord. This difference In the magnitudes of displacement thick- 
nesses due to Ice affects the pressure distribution and leads to a higher 
drag. 
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Fig. 25. Velocity profiles on the 
upper surface of Iced 
NACA 0012 airfoil, a = 
0.15° and R c = 1.5 x 10 6 . 



Fig. 26. Displacement thickness dis- 
tributions for clean and 
Iced NACA 0012 airfoil 
at a = -0.15° and R c = 1.5 
x 10&. 


To provide further Insight Into the behavior of the solutions, calculations 
were performed In which Invlscld and viscous flow equations were solved 
successively. That Is, rather than making several sweeps along the airfoil 
for a given pressure distribution, one boundary-layer sweep was made and 
the Invlscld flow was updated with the blowing velocity v n computed by 
the boundary- layer method. The results In Fig. 27 correspond to the vari- 
ation of the external velocity distribution on both surfaces of the airfoil 
with each cycle and show that the separation bubble of the previous calcu- 
lations becomes smaller and almost equal to that on the upper surface when 
the Initial Invlscld solution Is updated. In both cases, the separation 
and reattachment locations of the bubbles remain essentially unchanged 
after four cycles. As expected, the external velocity In the separated 
region Is uniform and decreases sharply near the reattachment point In 
accord with the behavior of separating and reattaching flows. The separa- 
tion bubble Is roughly ten percent of the chord, which Is In agreement with 
the experimental result of Bragg and Colrier. The lift Is nearly Independ- 
ent of the viscous effects for this angle of attack but the total drag 
coefficient requires the solution of the boundary- layer equations and the 
result Is different from that for a clean airfoil. 


Figures 28 to 30 show additional results for the same Iced airfoil. Com- 
parison of calculated lift and drag coefficients for a range of angles 
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Fig. 27. Variation of external velocity distribution with each cycle for 
Iced NACA 0012 airfoil at a = -0.15°, R c = 1.5 x 10&. 
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0.008 r- 



Flg. 28. Comparison of calculated and experimental results for the Iced 
NACA 0012 airfoil . R c = 1 .5 x 10 6 . 

of attack up to stall (a = 6°) Indicate that In general there Is good 
agreement with data. As expected, however, the behavior of the computed 
lift coefficients need Improvement. While the agreement at lower angles 
of attack Is satisfactory, It deteriorates with Increasing angle of attack 
due to the neglect of the wake effect. Figures 29 and 30 show the results 
for a - 6°. As can be seen from Fig. 29, at this angle of attack, there 
are substantial differences between the Invlscld and viscous velocity d 1 s - 
trlbutlon. Perhaps the most remarkable aspect of the calculations Is the 
behavior of the viscous flow solutions on the upper surface of the airfoil 
shown In Fig. 30a. The calculations Indicate approximately twenty-percent- 
chord leading-edge separation followed by fifteen percent marginally 
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Fig. 29. Comparison of Invlscld and viscous velocity distributions for 
the upper surface of the NACA 0012 Iced airfoil at a = 6°. 



Fig. 30. Variation of the local skin-friction coefficient Cf on the (a) 
upper, and (b) lower surfaces of the clean and Iced NACA 0012 
alrfol 1 at a = 6° . 


attached flow and sixty-five percent separated flow up to the trailing 
edge. It Is possible that with the Inclusion of wake effects, the extent 
of the flow separation on the airfoil will decrease. Nevertheless, from a 
numerical point of view, the calculations are able to cope well with such 
rather extensive flow separation. 
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